%% This file will produce fitted and countefactual expected win probabilities
%   for each project and aggregate them into expected danish and german maket shares
%    by project.
%  



%1. Read in theta from file (assume this has been estimated previously
%Contains 
% - sols, (each column is an estimate, presumably from a multi-start)
% - m, structure that contains model parameters
% - german_data, danish_data, so we can compute counterfactual distances

%% Either run these scripts or load their output files
%setup; 
%load data_structures;
%Main_mpec
%load est_point;

function [rho_nF, rho_nB, flag] = counterfactualSolve(m_nFC, thetahat, saveme) 

%% Counterfactual eliminating fixed costs of exporting to foreign market



%Now set up Jacobian pattern for the countefactual structure...
%These are the sparsity templates for each project the first row, is the 
%fringe share and the final column is theadding up constraint
J_m_ger = [[ones(1,m_nFC.num_firms_ger);eye(m_nFC.num_firms_ger)],ones(m_nFC.num_firms_ger+1,1)];
J_m_dnk = [[ones(1,m_nFC.num_firms_dnk);eye(m_nFC.num_firms_dnk)],ones(m_nFC.num_firms_dnk+1,1)];

%This is the Jacobian of the Constraints by the market share parameters
%(rhos)
JJ = [kron(eye(m_nFC.num_pro_ger),J_m_ger), zeros((m_nFC.num_firms_ger+1)*m_nFC.num_pro_ger, (m_nFC.num_firms_dnk+1)*m_nFC.num_pro_dnk) ; %These are the german market constraints
    zeros((m_nFC.num_firms_dnk+1)*m_nFC.num_pro_dnk, (m_nFC.num_firms_ger+1)*m_nFC.num_pro_ger), kron(eye(m_nFC.num_pro_dnk), J_m_dnk) ]; 



rho_ger = ones(m_nFC.num_pro_ger * (m_nFC.num_firms_ger + 1) ,1) ./ (m_nFC.num_firms_ger + 1);    
rho_dnk = ones(m_nFC.num_pro_dnk * (m_nFC.num_firms_dnk + 1) ,1) ./ (m_nFC.num_firms_dnk + 1);    
r0 = [rho_ger; rho_dnk];
Jac_Pattern_rho = JJ';
LB = [zeros(size(rho_ger,1)+size(rho_dnk,1),1)];
UB = [ones(size(rho_ger,1)+size(rho_dnk,1),1)];
ktropts = optimset('DerivativeCheck','off','Display','off',...
           'GradConstr','on','GradObj','on','TolCon',1E-6,'TolFun',1E-6,'TolX',1E-6,'JacobPattern',Jac_Pattern_rho);
[rho_nF, FVAL_nF, EXITFLAG_nF, OUTPUT_nF] = ktrlink(@(x_0) dummy_objective(x_0), r0, [],[],[],[],LB,UB,@(x_0) model_constraints(x_0, thetahat,m_nFC),ktropts,'knitro.opt');
% %[rho, FVAL, EXITFLAG, OUTPUT] = fmincon(@(x_0) dummy_objective(x_0), r0, [],[],[],[],LB,UB,@(x_0) model_constraints(x_0, thetahat,m),ktropts);

%% Counterfactual eliminating all boder-related costs of exporting

%Set the international variable border to the state international border...
thetahat_nB = thetahat;
thetahat_nB(end-1) = thetahat_nB(end);
[rho_nB, FVAL_nB, EXITFLAG_nB, OUTPUT_nB] = ktrlink(@(x_0) dummy_objective(x_0), rho_nF, [],[],[],[],LB,UB,@(x_0) model_constraints(x_0, thetahat_nB,m_nFC),ktropts,'knitro.opt');


flag = EXITFLAG_nF + EXITFLAG_nF;

%
%%%Optional save probably don't want it for the bootstrap...
if saveme
   save counterFact_out rho_nF EXITFLAG_nF OUTPUT_nF rho_nB EXITFLAG_nB OUTPUT_nB;
end